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ABSTRACT 


In this thesis Region of Interest Image Reconstruction from Projections 
using Wavelet transform has been studied. The convolution backprojection 
operator, used for reconstruction of images from projections, has a non-local 
filter which requires global projection data. This increases the objects’ ex- 
posure time to harmful radiation. It is proved that this filtering operation 
can be done on the chosen wavelet function instead of the projections. This 
localizes the filter and leads to the region of interest image reconstruction. In 
this thesis one-dimensional multiresolution analysis is used for region of in- 
terest image reconstruction and verified through simulations using computer 
generated images. The performance of the algorithm is compared with the 
other reported algorithms. 



Contents 


1 Introduction 1 

1.1 Review 1 

1.2 Organization of the Thesis 2 

2 Mathematical Preliminaries 3 

2.1 Image Reconstruction from Projections 3 

2.1.1 Radon Transform And Some Of Its Properties 3 

2.1.2 Inverse Radon Transform 5 

2.2 Wavelet Transform 8 

2.2.1 Continuous Wavelet Transform 8 

2.2.2 Discrete Parameter Wavelet transform 9 

2.2.3 Discrete Time Wavelet Transform and Discrete Wavelet 

Transform 11 

2.2.4 DPWT and multiresolution Analysis 11 

3 ROI Reconstruction Using Wavelets 15 

3.1 Wavelet Localization of The Radon Transform 15 

3.1.1 Nonlocality Of |if|W(0 15 

3.1.2 Wavelet Localization Of |^|W(^) 16 

3.2 Localization Of 1^1 W((^) With 2D MRA Reconstruction ... 17 

3.2.1 Multiresolution Reconstruction With 2D MRA 18 

3.2.2 ROI Reconstruction With 2D MRA 19 

3.3 Mutiresolution Reconstruction Using ID Wavelets 21 

3.4 ID MRA in Multiresolution Reconstruction 22 

3.4.1 Single Level MRA With CBP 22 



CONTENTS i, 


3 4.2 Determination of 2D filters 23 

3.4.3 Properties of the Constructed 2D filters 24 

3.4.4 Higher level Decomposition 25 

3.4.5 Region of Interest Reconstruction 26 

4 Implementation and Discussions 30 

4.1 Simulation Results 30 

4.2 Discussions 32 

4.3 Comparison of the Algorithms 32 

4.3.1 Aliasing effect 33 

4.3.2 Computational Complexity 35 

53 


5 Conclusion 



List of Figures 


2.1 Schematic for Radon transform operation 4 

2.2 Various filters, (a) Ram-Lak. (b) Shepp-Logan. (c) 

Lowpass cosine, (d) Generalized Hamming 7 

2.3 Computing wavelet coefficients through subband technique . . 13 

2.4 Inverse wavelet transform through subband techniques .... 13 

3.1 A simple combination of CBP and MRA 18 

3.2 Modification of the simple ‘CBP-MRA’ combination by mov- 
ing CBP inside 19 

3.3 Tiling of the frequency plane by the 2D wavelet transform, 

and the central slice used in the modified ID filter 20 

3.4 A simple combination of ID MRA and CBP 22 

3.5 Modification of the simple ‘MRA-CBP’ combination by mov- 
ing CBP inside 23 

3.6 Tiling of the frequency plane by the constructed 2D wavelet 

transform 24 

3.7 An equivalent implementation of MRA 25 

3.8 Tiling of the frequency plane by the constructed 2D wavelet 

transform 26 

3.9 Multiresolution reconstruction for three level MRA 27 

3.10 Region of reconstruction for three level MRA 29 

4.1 Sampling scheme for the region of interest reconstruction inputs. 31 

4.2 The aliasing effect produced in 2D MRA reconstruction for a 

particular angle 6 33 



LIST OF FIGURES 


IV 


4.3 The aliasing effect produced in 2D MRA reconstruction shown 

in two-dimensional frequency plane 34 

4.4 Shepp-Logan Phanthom head image used for simulations. ... 36 

4.5 Frequency responses of the ID analysis filters 36 

4.6 Frequency response of the 2D synthesis filter A(ni,n 2 ) con- 
structed with misaligned scaling function 37 

4.7 Frequency response of the 2D synthesis filter g{ni,n 2 ) 38 

4.8 Frequency response of the 2D synthesis filter hg{ni,n 2 ) 39 

4.9 Frequency response of the 2D synthesis filter /i/i 5 '(ni,n 2 ). . . 40 

4.10 Frequency response of the 2D synthesis filter /iAA(ni, 722 ). ... 41 

4.11 Illustration of the effect of misalignment of projection origin. . 42 

4.12 Multiresolution reconstruction: Image at resolution 2“^. ... 43 

4.13 Multiresolution reconstruction: Image at resolution 2“^ 43 

4.14 Multiresolution reconstruction: Image at resolution 2“^ 44 

4.15 Multiresolution reconstruction: Coarse image at resolution 2”^. 44 

4.16 Multiresolution reconstruction: Overall added image 44 

4.17 Frequency response of the 2D synthesis filter gR{ni,n 2 ) 45 

4.18 Frequency response of the 2D synthesis filter hyi{(ni, 712 ). ... 46 

4.19 Frequency response of the 2D synthesis filter 712 ). . . 47 

4.20 Centered region of interest image at resolution 2“^ 48 

4.21 Centered region of interest image at resolution 2~^ 48 

4.22 Centered region of interest image at resolution 2"^ 49 

4.23 Region of interest reconstruction: Coarse image at resolution 

2-3 49 

4.24 Overall added centered region of interest image 49 

4.25 Offcentered region of interest image at resolution 2“^ 50 

4.26 Offcentered region of interest image at resolution 2“^ 50 

4.27 Offcentered region of interest image at resolution 2-3 51 

4.28 Overall added offcentered region of interest image 51 

4.29 Computational requirements in number of multiphcations for 

the varying region of interest samples 52 



Chapter 1 


Introduction 


1.1 Review 

One of the conventional techniques used for the two-dimensional image recon- 
struction from its projections is the method of convolution back projection. 
In this method projections are filtered with a non-local filter which requires 
global projection data. Recently a new approach using wavelets has been 
reported in [9, 10], [13]-[18]. In reference [9], a wavelet localization of Radon 
transform is carried out and it is shown that wavelets with several zero mo- 
ments localizes the ramp filter used in convolution backprojection. This led 
a new approach for region of interest reconstruction of an image. 

Considering the drawbacks of the algorithm given in [9], Daleney et.al. 
proposed a new algorithm using the two-dimensional multiresolution analysis. 
The main goal of their algorithm is to reconstruct the image at various levels 
and a simple modification in the scheme results in the region of interest 
image reconstruction with reduced exposure. By introducing a set of angle 
dependent filters, derived from the 2D-MRA analysis filters, the convolution 
backprojection operation is shifted inside the 2D-MRA. The angle dependent 
filters are now combined with the non-local filters to obtain local filters. With 
this combined filters the region of interest is made possible. This is verified 
by simulation. 

In [11], the reconstruction of an image from one-dimensional wavelet 
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transform of the projections has been studied. It has been proved that in- 
verse Radon transform of wavelet transform of the projections at a particular 
scale is the two-dimensional wavelet transform of the given image. The two- 
dimensional wavelet is the inverse Radon transform of the one-dimensional 
wavelet operated on the projections. The algorithm was implemented in 
Fourier domain. 

In the present thesis work one-dimensional multiresolution analysis has 
been used for image reconstruction at various resolution levels and for region 
of interest image reconstruction. The underlying theory has been worked 
out in detail and explained with illustrative figures. The computer generated 
Shepp-Logan phanthom head has been used for implementing the algorithm 
and the results are discussed. The limitations of the proposed algorithm has 
been analyzed using the results obtained. 


1.2 Organization of the Thesis 

The thesis is organized as follows. 

• In chapter 2, the relevant results from Radon transform and Wavelet 
transform are reviewed. 

• In chapter 3, the reported works axe reviewed briefly and the proposed 
algorithm is explained in detail. 

• In chapter 4, the practical difficulties in implementing the algorithm is 
discussed. The simulation results for the computer generated images 
are shown. A discussion of the limitations of the algorithm is given. A 
comparison of our algorithm with the previously reported algorithms 
are given. 

• In chapter 5, we conclude the thesis with summarizing the algorithm’s 
performance. 



Chapter 2 


Mathematical Preliminaries 


This chapter briefly reviews the mathematical preliminaries used in forth- 
coming chapters. For detailed treatment about Radon transform readers are 
referred to [2] and [1], and for wavelets to references [5], [3], [6], and [4]. 


2.1 Image Reconstruction from Projections 

Image reconstructions from projections is the process of producing an im- 
age of a two-dimensional distribution (of some physical property) from the 
estimates of its line integral along a flnite number of hues at known loca- 
tions. The basic mathematical framework common to a large class of such 
reconstruction problem was first developed by Johann Radon, and is re- 
ferred to as Radon transform. Radon transform is a process which maps an 
n-dimensional function in Euclidean space 'RP' to a set of n — 1-dimensionaJ 
functions in Euclidean space along all possible angles. Here we restrict 

ourselves to Radon transform and its inverse for two-dimensional functions 
only. 

2.1.1 Radon Transform And Some Of Its Properties 

Definition : Let (x,y) designate coordinates of a point in a plane. Consider 
an arbitrary function f{x,y) defined over some domain D in R?. Then its 
Radon transform (also known as the projections), p{s,6) is given by 
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/ oo yoo 

/ y)S{xcos 6 + ysin 6 — s) dx dy (2.1) 

-OO J —oo 


If / is continuous and has compact support, then p(s, 9) uniquely determines 

/(a;, 2/) [!]• 

The function p(s, 9) at any particular value s and 9 is the line integral 
along the line s = x cos 9 y sin 9. In eqn. (2.1) d(-) is Dirac delta function, 
defined as 

1 if a: = 0 , . 

S(x) = < (2-2) 

0 if X 0 

Fig. 2.1 illustrates the way a two-dimensional function or image f{x,y) is 
mapped to p( 5 , 6 ). Practically 0^ is found only for finite number of angles 
and finite number of samples in region of support. 

Some properties of the Radon transform which we will refer to in this work 
are given below. In the following properties, let and p 2 {sy 6 ) 

represent the Radon transform of the functions fi^x^y) and / 2 (^??/) 
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respectively. is the Fourier transform of p{s,0) and .F(6,6) is the 

two-dimensional Fourier transform of f{x,y). 

1. Projection Slice theorem : The one-dimensional Fourier transform on 
s of p(s,d) is equal to the central slice, at an angle 9 of the two- 
dimensional Fourier transform of the object f{x,y). 

p{(,9) = cos 9,^ sm9) (2.3) 

2. Radon transform is linear, i.e., given the scalars, ai,a 2 G TZ, we have, 

R. X 

<^ifiix,y) + a2f2ix,y) < — > aipi(s,9) + a2P2{s,9) 

3. Convolution property : The Radon transform of the two-dimensional 
convolution of two two-dimensional signals fi{x,y) and f 2 ix,y) is the 
one-dimensional convolution of the Radon transform of the two signals. 

fi(x,y) (g)2 f 2 {x,y) M Pi{s,9) P2{s,9) (2.4) 

where (8)i represents one- dimensional convolution on ‘s’, ®2 represents 
two-dimensional convolution and R.T represents the Radon transform 
operation. 

2.1.2 Inverse Radon Transform 

Several inverse operations for obtaining the two-dimensional function from 
its projections, are reported in literature (refer [2] for them). Among them 
we consider only a particular method known as convolution back projection. 
Connected with that is the Back projection operator which operates on the 
projection p{s,9) to give a function b(x,y) as 

b(x,y) = ( p{xcos9 + ysm9.,9)d9 (2-5) 

Jo 

Back projection represents the accumulation of the ray sums of all the rays 
that pass through the point {x,y). Backprojection operator is not the inverse 
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of the Radon transform but the inverse Radon transform operation involves 
this operator. The inverse Radon transform is, 

1 /‘TT /*00 

fix, y) = ^ I di d6 (2.6) 

where p(^, 6) is the one- dimensional Fourier transform of p(s, d) with respect 
to ‘s’. We note that in equation (2.6), the inner integral represents the inverse 
Fourier transform of |^|p(i^, 0), i.e , it represents the filtering of the projections 
p(s, 9) with the ramp filter whose frequency response is |^|. The outer integral 
represents the backprojection operation on these filtered projection. 

In discrete domain implementation, |^| needs to be sampled to obtain 
a discrete filter. Since |^| is of infinite bandwidth and the value at higher 
frequency is much larger than at lower frequency, the aliasing effect will be 
enormous. This leads to bandlimiting |^| before sampling it. The window 
used IS rectangular window called the Ram-Lak filter. 

The Ram-Lak filter amplifies the high frequency components of the pro- 
jections. Since the noise strength in the high frequency region is generally 
more than the signal strength, we use a window with a gentle roll off char- 
acteristics to avoid the ringing at the edges. The characteristics of some 
functions that are normally used along with |^| are given below 


1. Ram-Lak 


Frequency response 


Impulse response 


Discrete impulse response 


Shepp-Logan 



l^irect(^d) 

^^[2sinc(2^o5) - sinc^(^O'S)] 

m = 0 
m 0 


4d’ 


^(7rm/2) 


K^ra^d 


Frequency response Hsl{0 

Impulse response ^sl(5) 

Discrete impulse response hsi^im) 


l^lsinc(<^d)rect(<^d) 

2( l+sinC27r$0S ) 
7r2(d2-4s2) 

2 

7r£i(l— 4m^) 
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Figure 2.2: Various filters, (a) Ram-Lak. (b) Shepp-Logan. (c) 

Lowpass cosine, (d) Generalized Hamming. 


3. Low-pass cosine 

Frequency response = |^| cos(7r^d)rect(^(i) 

Impulse response — f ) + WC'S + f )] 

Discrete impulse response hyf{m) = | — 2 ) + + 1)] 

4. Generalized Hamming for an o: : 0 < a < 1 

Frequency response ■^h(0 “ 

1^1 [cr -)- (1 — a) cos(27r^d)] rect(^d) 
Impulse response = 

oiKh{s) + ^ [/irlC-s -d) + ^h.l(s + d)] 

Discrete impulse response h-^{rn) = 

OiKhifn) + ^ [KT.{m - 1) -P K^{m + 1)] 

An additional problem associated with the is the need for filters 

with large lengths in spatial domain. This is so because each of these filters is 
non-differenciable at the origin of the frequency domain. The consequences 
of this and related effects are discussed in chapter 3. 
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2.2 Wavelet Transform 

Wavelets are relatively recent development in applied mathematics. The 
concept can be viewed as a synthesis of ideas which originated in diverse 
areas such as, physics, engineering, and pure mathematics. As a result, 
wavelets appeal to engineers as well as to scientists. Image coding, speech 
processing, coherent state study in physics, are but few of the applications 
of wavelets. Since there are four different types of wavelet transforms and as 
there are yet no standard names for them, the terminology in the literature 
can at times be confusing. We adapt the terminology used in [5]. 

2.2.1 Continuous Wavelet Transform 

The Fourier transform of a given signal gives the frequency content of the 
signal. However it fails to provide information about the time evolution of 
the frequency content, i.e., how the frequency contents of the given signal 
change with time. This time-frequency localization is achieved by the wavelet 
transform. The wavelet transform of signal f{t) evolving in time depends on 
two variables, scale ‘a’ and time ‘fe’. 

1 ft — h\ 

(CWT*/)(«, i) = IJjvJ j_^ fm j (2.7) 

where tp{t) is the mother wavelet function which satisfies the admissibility 
condition given by 

= 27r < oo (2.8) 

J-oo 1^1 

where is the Fourier transform of This requires 

roo 

/ dt = 0 

Joo 

and 

1^(01 < oo- 

When the wavelet function satisfies above condition, the inverse wavelet 
transform exists and is given by 
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The scaling used in the wavelet transform zooms into the short lived high 
frequency signals or zooms out the lengthy low frequency signal, by varying 
width of the support of the function This property of the wavelet 

transform distinguishes it from other time frequency analysis techniques such 
as the windowed Fourier Transform or STFT. 

2.2.2 Discrete Parameter Wavelet transform 

In eqn. (2.10), both ‘a’ and ‘6’ are continuous variables and there is a re- 
dundancy in the representation. The function f{t) can be recovered from 
{W^f){a, 6); for discrete values of a, b, if they are sampled sufficiently dense 
to satisfy the Nyquist criterion. As most of the common signals have short 
lived high frequency components and lengthy low frequency components, 
sampling the variables ‘a’ and ‘6’ will be optimal if ‘a’ is sampled in log scale 
and ‘6’ is made proportional to ‘a’ [4], i.e. a = a'^ and b = nboo^. So the 
discrete parameter wavelet transform is defined as 

(DPWT^/)(m, n) = -^ r /(t)V> (a^i - nbo) dt (2.10) 

ao' 

This can also be written as 

(DPWT^/)(m,n) = j (2.11) 

where It is of interest to study different wavelet 

functions and the constants oq and bo that permit perfect reconstruction 
(see the discussion in next paragraph also) of the function f{t). That is, 

fit) = EPPWT^/)(m, n)V>^„(t). (2.12) 

m,n 

Suppose there is over sampling on a and b, then eqn. (2.12) will not only 
hold but moreover, it is possible to have a non-unique representation of f{t), 
with respect to the same ^/’(t) [5]. 

The oversampling of the variables a and b leads to a new approach, rep- 
resentation of f{t) with non-orthonormal functions. The theory of frames 
in vector space is a generalization of the orthonormal decomposition of 
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functions, which includes non-orthonormal decompositions also. In an n- 
dimensional vector space V, a set of vectors {uj; z = 1 to m; m < n is 
considered. Any vector x G V is decomposed by {u,} to the inner product of 
X with every u,, represented as {x,v,). The set {uj} forms a frame in vector 
space V if and only if, for any non-zero vector x e V, 

m 

< B\\xf (2.13) 

2=1 

where A and B are constants, dependent on {uj only, called frame bounds, 
with 0 < A < 5 < oo. For a frame {u,}, a dual frame {u,} can be found so 
that any vector x E V can be obtained back exactly from its decomposition 
by {uj, as given by, 

m 

X = Yl{x,v^)v^. (2.14) 

2 = 1 

Some basic results from the theory of frames are listed below. 

1. A frame is tight \i A = B. 

2. If removal of a single element in a tight frame violates the lower bound 
of eqn. 2.13 then the tight frame is an exact frame. 

3. Elements of a tight, exact frame with A = B = 1 form an orthonormal 
basis. 

Eqn. (2.12) is a special case and another possible perfect reconstruction 
of the signal from eqn. (2.11) is also available. If are elements of a 

frame then a dual frame 'tpmn{t) can be found such that 

fit) = X^(DPWT^/)(m, n)iPr^4t). (2 15) 

m,n 

If ifrrin (t) forms a frame, it obeys the inequality 

-4||/Wf < EEK/W.V'mUOlP < -BII/WII" (2-16) 

m n 

with 0 < A < jB < oo. The constants A and B are frame bounds and are 
dependent only on ipmnit), and 

11/(011^ = / \fit)?dt < oo. 



2.2 Wavelet Transform 


11 


For some very special choice oi ip, and ho, it has been proved that if we 
choose ao = 2 and bo = 1 then {ipmn ’■ rn,n £ Z} constitutes an orthonor- 
mal basis of L^(TZ) ([3]). Amongst these functions, Daubechies orthonormal 
wavelets are well known and often used. In this thesis we use Daubechies 
orthonormal wavelets for simulation purpose. 

2.2.3 Discrete Time Wavelet Transform and Discrete 
Wavelet Transform 

For discrete time signals two different types of wavelet transforms are defined. 
One is the discrete time wavelet transform defined as, 

(DTWT/)(m, n) = ao ^ - nbo). (2.17) 

k 

This is time discretization of the discrete parameter wavelet transform. The 
other transform is the discrete wavelet transform defined as, 

(DWT/)(m,n) = 2-/2;^ miP{2-^k - n), (2.18) 

k 

where the discrete wavelet function ip{k) can be, but need not necessarily, a 
sampled version of its continuous counterpart. For detailed analysis of these 
two transforms readers are referred to [3, 5]. 

2.2.4 DPWT and multiresolution Analysis 

The idea of multiresolution analysis is very similar to dyadic subband de- 
composition and coding. In that, a lowpass filter h{n) and a highpass filter 
g{n) operates seperatly on the signal. After decimation by 2, the signal is 
decomposed into subband components. The same scheme is repeated on the 
lowpass output of the previous level to obtain the next level decomposition. 
It is well known that with proper choice of h(n) and giji) filters the signal 
can be perfectly reconstructed with upsampling and hin) and g{n) at every 
level [12]. This is related to discrete parameter wavelet transform decompo- 
sition. The theory developed by Mallat [6] is briefly reviewed below. 
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A sequence of vector spaces {V^ ■. j e Z) \s said to be a Multiresolution 
Analysis of if they satisfy the following properties and the theorems 

stated below [6]. 

i. . . . C F2 C Vi C Vo C Vli c y_2 C . . . 

ii. = L\1l) 
hi. ^3 = {0} 

iv. if fit) € Vj then /(2^t) € y, 

V. if fit) 6 y then /(t — A:) G y for ^ G 2 

Theorem : 1 For a given multiresolution analysis of L^iTl), there exists 
a function fit) such that {fmnit) ■ n E Z} forms a basis of Vm, where 
fmnii) = — n). The function fit) is unique and is referred to 

as a scaling function. 

Theorem : 2 Associated with each MRA, is a unique orthonormal wavelet 
function fit). 

Any signal /(t) G is projected onto the vector space Vq by 

dt (2.19) 

A=nT 

This^(n) is decomposed using filter gin) to get discrete parameter wavelet 
transform of fit) at level. The decomposition by ^(n) results in the 
projection of fit) onto y, denoted as ci. Further decomposition of ci with 
gin) and h(n) results in c„ and discrete parameter wavelet transform of fit) 
at 2”” level. 

The filters gin) and ^(n) are related to wavelet and scaling functions and 
given by, 

hin) = ^ / fit)fi‘^'t - n) dt 
gin) = \j fii)fi‘^^ - «) dt, 



( 2 . 20 ) 

( 2 . 21 ) 
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Figure 2.3: Computing wavelet coefficients through subband technique 



Figure 2.4: Inverse wavelet transform through subband techniques 


and 


p-i 


4>{t) = 

1=0 

(2.22) 

m = 

1=0 

(2.23) 

From the above theory it is noted that 


4(n) = J fit)i^kni't)dt. = DP'WT^{k,n) 

(2.24) 

Eqn. (2.24), gives a simple algorithm to obtain discrete parameter wavelet 

transform of the signal. The decomposition scheme is shown 
reconstruction scheme is given in Fig. 2.4. 

in Fig. 2.3. The 
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In practice discrete time signals are needed to be analyzed rather than 
continuous time signals. So the discrete time signal is assumed as f{n). As 
long as we can reproduce f{n), the corresponding f{t) is not coming m to the 
picture. In this thesis we use this multiresolution analysis for discrete param- 
eter wavelet transform using Daubechies-10 compactly supported wavelet. 



Chapter 3 


ROI Reconstruction Using 
Wavelets 


In this chapter we consider the region of interest reconstruction from mini- 
mum possible projections. This problem has been considered and tried out by 
few authors also in different ways [13]-[18],[10, 9]. Among them we consider 
one particular method[9], which combines the non-local filter with wavelet 
function to obtain a local filter. In section 3.1 and 3.2 we review the methods 
from references [9] and [10] respectively. Section 3.3 reviews the theory of 
reference [11]. Section 3.4 discusses the proposed algorithm. 


3.1 Wavelet Localization of The Radon Trans- 
form 

In this section, the reason for the non-locality of the filter is first 

analysed. The method to overcome this problem, as given in [9], is briefly 
discussed. 

3.1.1 Nonlocality Of |^|1T(^) 

As mentioned in the chapter 2. the filter is not locally supported. 

This stems from the fact that the function |^| is not differentiable at the 
origin. The non-locality implies that local calculation of the convolution 
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in eqn. (2.6) will require global values of the Radon transform. The non- 
differentiability at the origin can not be significantly altered by the window 
TT(^), without distorting the structure of the image. 

Consider the convolution, (which is multipUcation in the frequency do- 
main) as a filter operation denoted by and we will write it as a two step 
process. The first step is multiplication by in the frequency domain (or 
differentiation in the time domain). The second step is the application of 
Hilbert transform H. Thus the filtering process is given by 

= IFT{isign(nK)/(f)} 

= IFT{|{|/({)}. (3.1) 

where 

1 if re > 0 ^ ^ 

sign(0 = { ^ (3-2) 

0 if a: < 0. 

The differentiation operation is local, but the Hilbert transform is not, since 
it imposes non-differentiability upon the Fourier transform of any smooth 
function whose average value is not zero. 

3.1.2 Wavelet Localization Of 

The imposition of non-differentiablity at the origin in the frequency domain 
and the subsequent spreading of the support of the function, will not occur 
if the function’s Fourier transform is zero at the origin and the function 
has several zero moments. But we can not impose conditions on the Radon 
transform of an image. The trick here is to combine the filter with 

some other function which has the desired properties as well as provides 
a method of reconstruction with the combined filter. This can be done as 
follows. 

Consider the convolution backprojection equation (2.6) given by, 

= /?(»,«) 


(3.3) 

(3.4) 
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m^n 

(3.5) 

I Y^Cm,n{0)^{lpm,n{s)}de 

m,n 

(3.6) 


This gives a novel reconstruction formula where the filtering is done on 
wavelet function. By choosing a wavelet function with several zero moments 
we can localize the and hence the inverse Radon transform. When 

1 pm, nis) is compactly supported, Cm,n{0) can be calculated using local infor- 
mation. In this algorithm to obtain the region of interest reconstruction 
the projections are taken along the region of interest with a small offset 
accounting for both the calculation of Cm,n(^) and the inverse operation. Im- 
plementation of this is discussed in [9]. The major drawbacks [9, 10] in their 
algorithm are listed below. 

1. Though, the reconstruction at the region of interest is good, the recon- 
structed image has severe artifacts outside the region of interest. 

2. The algorithm is prone to aliasing error when the region of interest is off 
centered. The process of centering the projections before interpolation 
and reconstruction introduces this aliasing error. 

3. It is not mentioned in the reference, how multiple non-concentric re- 
gions can be handled. 

3.2 Localization Of |e |w (^) With 2D MRA 
Reconstruction 

In [10], Delaney and Bresler uses the basic ideas from [9] and taking into con- 
sideration the algorithm’s drawbacks, proposes a new algorithm which com- 
bines the filter with two-dimensional separable MRA analysis filters. Their 
algorithm is briefly discussed. 
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Figure 3.1: A simple combination of CBP and MRA. 


3.2.1 Multiresolution Reconstruction With 2D MRA 

Consider Fig. 3.1. Since the MRA gives perfect reconstruction f{ni,n 2 ) and 
the synthesised output will be exactly equal. The analysis is carried out in 
two-dimension as 

/i(ni,n2) = /(rai,n2) ®2d hh{ni,n2) (3.7) 

By convolution property of Radon transform, the equation is modified to 

fi{ni,n 2 ) = IRT{p(n, 9) 0id khkin)} (3.8) 


where, 

hhiO, n) = IFTjdiM^ cos e)h{^ sin 0)} (3.9) 

This modified scheme is shown in Fig. 3.2, in which ‘/F blocks represents 
the one-dimensional angle dependent filters obtained from eqn. (3.9). The 
advantage of having this representation is that the convolution part of CBP 
can be now combined with the ‘A:’ filters. As proved in [9], the combination 
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Figure 3.2: Modification of the simple ‘CBP-MRA’ combination by moving 
CBP inside. 

will result in local filters if the filters have zero at the origin and in the 
neighbourhood of origin. 

3.2.2 ROI Reconstruction With 2D MRA 

By examining how the two-dimensional separable wavelet transform tiles the 
frequency plane, we can get a clear perception of the properties of the one- 
dimensional angle dependent filters. The tiling is shown in simplified form 
in Fig. 3.3. Note that spectra of actual filters overlap with each other to 
some extent. By the projection slice theorem, the one-dimensional filter at 
particular angle has a passband from to So for all 6 the k^hiO, n) are 
lowpass filters. The reconstruction obtained from this particular channel is 
of very low bandwidth. This reduces the number of projection angles needed 
in Radon transform to a minimum. So projections taken at less number of 
angles are given to this channel, which means less exposure of the object to 
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Figure 3.3: Tiling of the frequency plane by the 2D wavelet transform, and 
the central slice used in the modified ID filter. 

radiation. 

On the contrary, other channels use high-pass filters of different nature 
each having zero at the origin of frequency plane and zero in the neighbour- 
hood. So if we combine the function |^ilT(^) with these filters they result 
in compactly supported filters as explained earlier. So with more angular 
sampling, the projections are taken only around region of interest with an 
offset accounting for the filter length. The over all algorithm is 

• Find the angle dependent filters using eqn. (3.8). 

• Convolve these filters with 

• Input the projections corresponding to a reduced number of angles to 
the lowpass channel. 

• Input the projections at all angles around the region of interest to other 
channels. 


Backproject the data on all channels. 
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• Downsampling and upsampling is carried out along row and column. 

• 2D MRA synthesis is carried, out in every channel. 

• The image at various resolution levels are found. Add the images at 
all levels, the result is a ROI reconstructed image. 


3.3 Mutiresolution Reconstruction Using ID 
Wavelets 


In [11] Peyrin et.al. studied the reconstruction of tomographic images from 
the one-dimensional wavelet transform of its projections. They achieved 
this in Fourier domain only. They proved that, the reconstructed image 
from the wavelet transform (at a fixed scale) of the the projections is the 
two-dimensional wavelet transform of the original image at that scale. A 
modified proof is also given in [12]. This result is important because it gives 
the procedure to construct a two-dimensional wavelet function from the one- 
dimensional wavelet function. We use this theory in the algorithm developed 
in this thesis. 

The two-dimensional wavelet function constructed by inverse Radon trans- 
forming the one-dimensional wavelet function is different from the conven- 
tional separable wavelet functions used for images. It is isotropic in frequency 
domain and non-separable. The implementation scheme given in [11] is re- 
produced here. 


• For an assumed one-dimensional wavelet function, a two-dimensional 
wavelet function is found using the Fourier transform method. 

• The Fourier transform of the projections are taken. 

• The projections axe multiplied with the wavelet function in Fourier 
domain to find the wavelet transforms of the projections. 


Using Fourier domain reconstruction the wavelet transformed projec- 
tions are inverse Radon transformed. 
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Figure 3.4: A simple combination of ID MRA and CBP. 


• The output is multiplied with the two-dimensional wavelet function in 
Fourier domain to find the images at various resolution levels. 

3.4 ID MRA in Multiresolution Reconstruc- 
tion 

In this section using the proof given in [11] a new algorithm is proposed 
and discussed. Unlike the work described by Delaney in [10], the proposed 
method primarily uses the ID-MRA analysis of the projections. At the syn- 
thesis end, the convolution backprojection operator is applied before the 
synthesis filters, in each channel respectively As an example a simple one 
level MRA decomposition is taken to prove the result. Then the scheme is 
extended for higher level of decompositions. 

3.4.1 Single Level MRA With CBP 

As shown in Fig. 2.3, the projection p{s,9) is first given as an input to 
a single level one-dimensional MRA filter bank. The decomposed signals, 
Ci{n,9),di{n,0) are used as inputs to the synthesis stage as shown in Fig. 
3.4. Since multiresolution analysis and synthesis gives perfect reconstruction 
the p{s, 9) will be exactly reconstructed at the synthesis end. Define a two- 
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Figure 3.5: Modification of the simple ‘MRA-CBP’ combination by moving 
CBP inside. 

dimensional function /i(n, 0) as 

h{n,6) — h{n) for every 6, (3.10) 

The synthesis equation can be written as 

P 2 {n, 6) = pi{n, 6) h{n, 9) (3.11) 

By convolution property of Radon trans form we have, 

/2(ni, raz) = /i(ni, nj) ®2d h(ni, 712 ) (3.12) 

Similarly, in other channels, the convolution backprojection operator is is also 
moved in before the synthesis filters. The modified implementation scheme 
is shown in Fig. 3.5. 

3.4.2 Determination of 2D filters 

Though the two-dimensional filters are inverse Radon transform of the cor- 
responding one-dimensional filters, the straight forward implementation of 
convolution backprojection on one- dimensional filters will result in filters 
having extremely poor frequency response, because of reduced number of 
samples. Here we propose an alternative scheme, described below. 

• From h{n) and g(n) find the corresponding one-dimensional wavelet 
and scaling functions using eqn. (2.24). This can be done at any re- 
quired sampling rate depending on the accuracy needed. 
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Figure 3 . 6 : Tiling of the frequency plane by the constructed 2 D wavelet 
transform. 


• Find the two-dimensional wavelet and scaling functions by inverse Radon 
transforming the one-dimensioneil wavelet and scaling functions respec- 
tively. 


• Find the two-dimensional filters using 

/i(ni,n2) = J J (f>2d{x,y)(l>2di2x -ni, 2 y -n2)dxdy (3.13) 
g(ni,n2) = j J i^2dix,y)<l>2d{2x - ni, 2 y - n2)dx dy (3.14) 


This is done only once in the entire process of first reconstruction. After that 
for repeated reconstruction, for the same wavelet function, we use the same 
filters. 


3.4.3 Properties of the Constructed 2D filters 

As mentioned before the properties of the two-dimensional filters affect the 
performance of region of interest reconstruction. Since h{ni^n<^ is inverse 
Radon transform of h{n) we can perceive it as an isotropic lowpass filter. 
Si mil arly g{n\,n2) has the passband as shown in Fig. 3 . 6 . We note that 
gipi, 712) is zero at the origin and in the neighbourhood in the Fourier domain. 
From the discussion given earlier, if |^|W (0 i® combined with this filter, it 
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Figure 3.7: An equivalent implementation of MRA. 

will result in a smaller length filter. This is the major advantage of this 
modified technique. 

3.4.4 Higher level Decomposition 

Consider the MRA synthesis filter bank shown in Fig. 2.4. The convolution 
backprojection operator cannot be shifted before the upsampling operations. 
So we face difficulties at higher level of decompositions. Let us view this 
problem in the analysis domain. Consider a three level analysis filter bank 
for the one-dimensional MRA. The filters which are present after the down- 
sampling, can be placed before the downsampling block, thereby resulting in 
upsampling this filter. The modified analysis filter bank is shown in Fig. 3.7. 
But in the two-dimensional case, as explained earlier, we need an explicit ex- 
pression to calculate two-dimensional filter from the two-dimensional wavelet 
and scaling functions. Using the theory of muliresolution analysis, the filters 
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Figure 3.8: Tiling of the frequency plane by the constructed 2D wavelet 
transform. 

can be obtained using the equations 

hk{n) = J - n)dt (3.15) 

h-igin) = J - n)dt (3.16) 

This can be extended to two-dimension as 

hkini,n2) = J J cf>2d{x,y)<f>2d{^’‘x-ni,2’'y -n2)dxdy (3.17) 
hk-ig{ni,n2) = j J 'tp 2 d{x, y) 4 > 2 d{. 2 ^x - ui, - n 2 )dx dy (3.18) 

The passband characteristics of the two-dimensional filters are shown in Fig. 
3.8. The convolution backprojection can now be carried out before the fil- 
tering operation. The modified implementation scheme for a three level mul- 
tiresolution reconstruction is shown in Fig. 3.9. 

3.4.5 Region of Interest Reconstruction 

Referring to Fig. 3.8 we observe that the filters hk-ig{ni,n 2 ) have zero value 
at the origin as well as in the neighborhood of the origin. So if we can 
combine 1^11F(0 with these functions in the two-dimensional case, we can 
obtain short length filter. The convolution backprojection operation can be 
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9di{n,6) CBP ► g(ni,n2) 

9d2{n,0) — CBP — ► hg{ni,n2) 

9dz{'<^i^) — CBP — ► hhg(ni,n2) 

9c:z{P'->^) — CBP — ^ hhh(ni,n2) 

Figure 3.9: Multiresolution reconstruction for three level MRA. 

represented as two blocks, first is the convolution with |i^|W(^) and next 
is backprojection operation. Using the convolution property of the Radon 
transform we prove here that the one-dimensional filter which is an integral 
part of convolution backprojection can be moved over the backprojection 
block, resulting in a two-dimensional filtering with the inverse Radon trans- 
form of the one- dimensional filter. Consider the convolution backprojection 
equation 

= j p{s,D)S:u 1 FTu{\(\W(()} de (3.19) 

= /lFT,,{p(e.*)}®iiIFT,4|f|iy(f)}i« (3.20) 

= / IFT,,|i|;|-|{!lF(f)| ®.^IFTw{|f|»^(«} d6(3.21) 

IFTi4|{HV(0}j ®wIPTij{|{|lV({)} de 
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Equation (3.23) is the inverse Radon transform of one-dimensional convolu- 
tion of two functions. Applying the convolution property of Radon transform 
it can be written as, 


fix,y) = 


/ 


IFT 


f 


02d (3.23) 


/ IFTi,{|eilF(0}®i<i IFTid{ieifF(0} d6. 


f(^,y) = 



®2J [v'S + g 


(3.24) 


where 

W 2 d{^x-, (y) = fF(^ = COS 0 4- iy sin 0). 

Using this approach we combine the |^|tU(^) in the two-dimension with the 
two-dimensional synthesis filters. This combination is advantaguous only 
when it is combined with highpass or bandpass filters. Further more it gives 
computational disadvantage with lowpass filter. Since the two-dimensional 
filter is non-separable, the filtering in two-dimension is computationally in- 
tensive than than filtering in one-dimension. Since the conversion of |^lfU(^) 
to two-dimension in one channel is independent of the other channels, this 
is not done in lowpass channel to avoid the computational cost. The Imple- 
mentation algorithm is clearly shown in the Fig. 3.10. 



3 4 ID MRA in Multiresolution Reconstruction 


29 


9di{s,9) 


gd2{s,9) 


9dn{s,0) 


5'cn('S) 9 ) 



Figure 3.10: Region of reconstruction for three level MRA. 









Chapter 4 


Implementation and 
Discussions 


The algorithm is implemented using Daubechies ten compactly supported 
orthonormal wavelets. The projections are simulated for the head phan- 
thom(Fig. 4.4). The results are given in the section 4.1. A brief discussion 
on the results is given in ??. A comparison of the new algorithm with other 
algorithms is given in section 4.3. 


4.1 Simulation Results 

The implementation scheme for multiresolution reconstruction is shown Fig. 
3.9 and the region of interest reconstruction in Fig. 3.10. 

As discussed, the one-dimensional scaling function and wavelet are in- 
verse Radon transformed to obtain the two-dimensional scaling function and 
wavelet. Using eqn. (3.18), /i(ni,n 2 ) and 5 '(ni,n 2 ) are obtained. To our sur- 
prise the filter responses are not what we expected or rather what is given 
in theory. This is due to the non-alignment of center of projections in the 
inverse Radon transform operation. To explain briefly, the center point of 
the projections should be aligned in backprojection, otherwise the function 
in two-dimension cannot be reconstructed exactly. 

To align the centers, zeropadding of the sequence is carried out and the 
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Figure 4.1: Sampling scheme for the region of interest reconstruction inputs. 

two-dimensional scaling function and wavelet are found. To illustrate this 
effect, for a misaligned scaling function, inverse Radon transform is taken 
to find the /i(ni,n 2 ). The frequency response of this filter is shown in Fig. 
??. To center them, we used 197 zeros padded on the left of scaling function 
and 16 zeros padded on the left of wavelet. This is found by repeated trials, 
to make the Fourier transform of the two-dimensional filters’ cross section 
correspond to one-dimensional filter response. The one-dimensional filters’ 
responses are shown in Fig. 4.5 and the two-dimensional filters responses are 
shown in Figs. 4.7 to 4.10. Since the convolution of projections with MRA 
filters introduce noninteger phase shift, the origin of the projections is shifted 
and the zeropadding procedure cannot solve this as the shift is a noninte- 
ger. The effect of this noninteger phase shift on the reconstructed image is 
illustrated for Shepp-Logan head phanthom using conventional convolution 
backprojection operator is shown in Fig. 4.11. To correct this problem, an 
offset value is added with the s = x cos ^ -f- J/ sin 9. The simulation results on 
Shepp-Logan head phanthom are shown in Fig. 4.20-4.24. 
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For region of interest reconstruction the scheme shown in Fig. 3.10 is 
implemented. The projections corresponding to a reduced number of angles 
are given input to the lowpass channel and the input to the other channels 
are illustrated clearly in Fig. 4.1. The region of interest reconstruction filter 
frequency responses are shown in Fig. 4.17-4.19. Simulation results are shown 
in Fig. 4.20-4.24. The simulation results, to demonstrate that the offcentered 
region of interest reconstruction is also possible, are shown in Fig. 4 25 to 
4.28. 


4.2 Discussions 

As we stated previously, the new algorithm for region of interest image re- 
construction from projections work fairly well at resolution levels and 
2~^. The edge information on these levels are fairly acceptable. The images 
at resolution levels 2"^ however is blurred as expected. The same is true for 
region of interest reconstructed images at various levels. 

The limitations of the algorithm are twofold. 

• The process of finding the two dimensional filters requires inverse Radon 
transforming the scaling and wavelet functions respectively. An align- 
ment of the centers for all angles has to be carried out before inverse 
Radon transforming. 

• As the one-dimensional multiresolution analysis filters introduce non- 
integer phase delays on the projections, an offset is added while taking 
the inverse Radon transform. 

4.3 Comparison of the Algorithms 

As discussed earlier, Olson’s algorithm, given in [9], is inferior in many as- 
pects as given in [10], we compare here our algorithm only with the results 
reported in [10]. 
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Figure 4.2: The aliasing elfect produced in 2D MRA reconstruction for a 
particular angle 0 . 

4.3.1 Aliasing effect 

The Delaney’s algorithm gives aliasing at high frequencies thereby degrading 
the image. This is due to the fact that the one-dimensional discrete frequency 
is restricted between — tt to tt. This can be explained from the frequency 
response plot of the two-dimensional MRA filters and its Radon transformed 
one-dimensional MRA filters. This is shown in the Fig. 4.3 and the aliasing 
effect for a single filter at an angle 6 is shown in Fig. 4.2. In lowpass filters no 
aliasing occurs. But in highpass filters when the projection angle is not equal 
to 0 or TT /2 the two-dimensional frequency plane extends from — tt/ cos 9 to 
tt/ cos (9 for 61 < f (and for 6 < f frequency plane is from -tt/ sin 61 to tt/ sin 9 . 
Since this is grater than the span from — tt to tt, the aliasing effect comes. 

As the images in general have very low power in the high frequency com- 
ponents, this was not observed in the results reported in [10]. In the proposed 
algorithm, as the synthesis filters are the inverse Radon transform of the one- 
dimensional filters, the synthesis filters have zero outside the circle of radius 
TT in two-dimensional frequency plane. This proves that the developed algo- 
rithm according to theory will give better results. 
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Figure 4.3: The aliasing effect produced in 2D MRA leconstruction shown 
in two-dimensional frequency plane. 
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4.3.2 Computational Complexity 

Designing the two-dimensional filters is computationally very intensive, as we 
need to find inverse Radon transform of the fine sampled scaling and wavelet 
functions. But this is to be done only once while setting up the algorithm. So 
this is not considered in the run time computational requirement calculations. 

The MRA synthesis is done in the two-dimensional domain in which the 
filters are non-separable. So the number of computations is increased to 

X Ni whereas in separable filter implementation it is iV x iVi (A^ is the 
length of the two-dimensional filter in both direction or one-dimensional filter 
length and Ni is sequence length). As mentioned before the lowpass channel 
is not modified in region of interest reconstruction, from the multiresolution 
reconstruction algorithm. This significantly reduces the number of multi- 
plications. The overall computational requirements for the multiresolution 
reconstruction is shown in Fig. 4.29. In that the number of samples con- 
sidered for projections is 100 and number of angles considered is 128. The 
image is sampled in 256 x 256 grid and filtered. The computational require- 
ments for region of interest reconstruction is shown in Fig. ??. The region 
of interest radius is 30 and the sampling scheme is shown in Fig. 4.1. The 
number of projection samples for lowpass filter is 256 and number of angles 
considered is 32. The image is sampled only in the region of interest. 

Note : In two-dimensional frequency response plots, the first plot denoted 
as (a), is with the hidden surfaces and the second one in the same figure, 
denoted as (b ), is without the hidden surfaces. 








Figure 4.7; Frequency resf 




Figure 4.8; Frequency response 


of the 2D synthesis filter hg{nx,n-i). 
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Figure 4.11: Illustration of the effect of misalignment of projection origin. 
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Figure 4.12: Multiresolution reconstruction- Image at resolution 2 


Figure 4.13: Multiresolution reconstruction: Image at resolution 2 ^ 
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Figure 4.14: Multiresolution reconstruction: Image at resolution 2 


Figure 4.15: Multiresolution reconstruction: Coarse image at resolution 2 ^ 


Figure 4.16: Multiresolution reconstruction: Overall added image. 
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Figure 4.17: Frequency response 


of the 2D synthesis filter gR{nr,n 2 ). 
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Figure 4.18: Frequency response 


of the 2D synthesis filter hgiiini,n 2 ). 
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Figure 4.20; Centered region of interest image at resolution 2 ^ 


Figure 4.21: Centered region of interest image at resolution 2 ^ 
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Figure 4.22: Centered region of interest image at resolution 2 


Figure 4.23: Region of interest reconstruction. Coarse image at resolution 

2 - 3 . 


Figure 4.24: Overall added centered region of interest image. 
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Figure 4.25: Offcentered region of interest image at resolution 2 ^ 


Figure 4.26: Offcentered region of interest image at resolution 2 ^ 



4 3 Comparison of the Algorithms 


51 


Figure 4.27. Offcentered region of interest image at resolution 2 


Figure 4.28: Overall added offcentered region of interest image. 
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Figure 4.29: Computational requirements in number of multiplications for 
the varying region of interest samples. 





Chapter 5 


Conclusion 


The objective of the thesis was to explore a new algorithm for region of in- 
terest image reconstruction. It has been found that the proposed algorit hm 
is able to extract the edge information at resolution levels 2 ~^ and 2 ~^. But 
the overall reconstruction is found to be visibly poorer than the conventional 
convolution backprojection method. This is due to the extensive operations 
done on the projection data and the inverse Radon transform of the signal 
which has 2” — 1 zeros at the level. The aim of the algorithm developed 
was not solely to reconstruct the image from projections, but also to reduce 
the amount of radiation exposure, i.e. to use minimum number of projection 
samples. By combining the non-local filter with the two-dimensional syn- 
thesis filters the localization is achieved. The algorithm is also applicable 
for offcentered region of interest image reconstruction.the present algorithm 
does not have the aliasing error unlike the algorithm given in [10]. This is 
achieved at the cost of increased computations. The algorithm is found to 
be suitable for reconstructing the region of interest images at 2“^ and 2“^ 
resolution levels. This is because the image reconstruction at these levels 
needs less computaion, at the same time gives the required edge information. 
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